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Hamilton-Jacobi  (H-J)  equations  are  frequently  encountered  in  applications,  e.g.  in  con¬ 
trol  theory  and  differential  games.  H-J  equations  are  closely  related  to  hyperbolic  conserva¬ 
tion  laws  f  m  one  space  dimension  the  former  is  simply  the  integrated  version  of  the  latter. 
Similarity  also  exists  for  the  multi-dimensional  case,  and  this  is  helpful  in  the  design  of  dif¬ 
ference  approximations.  In  this  paper  we  investigate  high  order  essentially  non-oscillatory 
(ENO)  schemes  for  H-J  equations,  which  yield  uniform  high  order  accuracy  in  smooth  re¬ 
gions  and  resolve  discontinuities  in  the  derivatives  sharply.  The  ENO  scheme  construction 
procedure  is  adapted  from  that  for  hyperbolic  conservation  laws.  We  numerically  test  the 
schemes  on  a  variety  of  one-dimensional  and  two-dimensional  problems,  including  a  problem 
related  to  control-optimization,  and  observe  high  order  accuracy  in  smooth  regions,  good 
resolution  of  discontinuities  in  the  derivatives,  and  convergence  to  viscosity  solutions. 
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also  supported  in  part  by  NSF  Grant  DMS  88-10150. 


1.  INTRODUCTION 


The  Hamilton- Jacobi  (H-J)  equation 

(11) 

l  <f>(x,0 )  =  4> 0(x) 

appears  often  in  applications,  e.g.,  in  control  theory  and  differential  games.  The  solutions 
to  (1.1)  typically  are  continuous  but  with  discontinuous  derivatives,  even  if  the  initial 
condition  4>o(x)  is  C°°.  The  non-uniqueness  of  such  solutions  to  (1.1)  also  necessitates  the 
introduction  of  the  notions  of  entrcoy  condition  and  viscosity  solutions,  to  single  out  a 
unique  practically  relevant  solution.  See,  e.g.  Crandall  and  Lions  [1]  for  details. 

An  important  class  of  numerical  methods  for  solving  (1.1)  is  the  class  of  monotone 
schemes  discussed  by  Crandall  and  Lions  [2].  Monotone  schemes  are  proven  convergent  to 
the  viscosity  solutions.  Unfortunately  monotone  schemes  are  at  most  first  order  accurate. 
Traditional  high  order  methods  are  unsuitable,  because  spurious  oscillations  will  generally 
occur  in  the  presence  of  discontinuous  derivatives. 


There  is  a  close  relation  between  (1.1)  and  a  hyperbolic  conservation  law 

|  ut  +  £?=  i  fi(u) x.  =  0 
l  u(x,0)  =  u0(x). 


(1.2) 


In  fact,  for  the  one  dimensional  case  d  =  1,  (1.1)  is  equivalent  to  (1.2)  if  we  take  u  =  <j>x. 
For  multi-dimensions  this  direct  correspondence  disappears,  but  in  some  sense  we  can  still 
think  about  (1.1)  as  (1.2)  “integrated  once”.  Hence  successful  numerical  methodology  for 
solving  hyperbolic  conservation  laws  (1.2)  should  also  be  applicable  to  the  H-J  equation 
(1.1). 
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Essentially  non-oscillating  (ENO)  schemes  have  been  very  successful  in  solving  the 
hyperbolic  conservation  law  (1.2),  Harten  and  Osher  [3],  Harten,  Engquist,  Osher  and 
Chakra varthy  [4],  Shu  and  Osher  [7],  [8].  The  key  idea  is  an  adaptive  stencil  interpolation 
•which  automatically  obtains  information  from  the  locally  smoothest  region,  hence  yields  a 
uniformly  high  order  essentially  non-oscillatory  approximation  for  piecewise  smooth  func¬ 
tions.  We  summarize  this  ENO  interpolation  procedure  as  follows: 


Given  point  values  f(xj),j  =  0, ±1,±2, ...  of  a  (usually  piecewise  smooth)  functions 
at  discrete  nodes  ij,  we  associate  a  r-th  degree  polynomial  P^y2(x)  with  each  interval 
constructed  inductively  as  follows: 


(1)  —  f[^j  “b  f[xj ,  Xj+\ ](^  ^min  — 3 


(1.3) 


(2)  If  in1J  and  pj£ \^{x)  are  both  defined,  then  let 

O  ^  —  /[x^i- 1) ,  •  •  •  ,  X .  (<->).*] 

min  min  ' 

b(f)  =  /[ ,... 

min  min  *  4 


and 


(i)  if  |«W|  >  \b(%  then  c<'>  =  b«\  _  1 

otherwise  c<'>  =  c('>,  *£??  = 


(») 


l('-D 


+/-1 


pf,t 
J+i/  2 


(*)  =  Pi 


/+  1/2(X)  + 


n<*> 


n  (x  -  z0 


i=k«-y 


□ 
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In  the  above  procedure  /[•, ...  ,•]  are  the  usual  Newton  divided  differences.  Notice  that 
we  can  also  start  from  one  node  Xj  to  build  a  polynomial  Pj'T{x)  using  the  same  procedure. 

In  [6]  Osher  and  Sethian  constructed  ENO  type  schemes  and  applied  them  to  a  class  of  H- 
J  equations  and  perturbations,  arising  in  front  propagation  problems.  They  achieved  very 
good  numerical  results.  In  this  paper  we  provide  a  more  general  ENO  scheme  construction 
procedure,  mainly  by  considering  different  multidimensional  monotone  building  blocks. 
We  then  numerically  test  the  schemes  on  a  variety  of  one  and  two  dimensional  problems 
including  a  problem  related  to  control-optimization,  check  the  accuracy  in  smooth  regions, 
the  resolution  of  discontinuities  in  derivatives,  and  the  phenomenon  of  convergence  to 
viscosity  solutions.  Concluding  remarks  are  given  in  Section  4. 


2.  SCHEME  CONSTRUCTION 


For  simplicity  of  exposition  we  take  d  =  2  in  (1.1),  and  use  x,  y  instead  of  Xj ,  x i'- 


<{>t  +  H{<t> Xi  <f>y)  =  0 

<f>(x,y,  0)  =  <f>0 (ar,y) 


(2.1) 


For  mesh  sizes  Ax,  Ay,  At,  <f>^  will  denote  a  numerical  approximation  to  the  viscosity 
solution  4>(xi,yj,tn )  =  <^(iAx,j  Ay,  nAt)  of  (2.1).  We  also  use  standard  notation  such  as 


\  A ±<f>ij  =  -  <t>ij ) 


We  start  with  a  first  order  monotone  scheme  [2]: 


^  =  4%  -  A tH( 


A-*r,  a  us 


\  Ax  ’  Ax  ’  Ay  ’  Ay 


) 


(2.2) 
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where  H  is  a  Lipschitz  continuous  monotone  flux  consistent  with  H : 


H(u,u,v,v )  =  H(u,v). 


Monotonicity  here  means  that  H  is  non-increasing  in  its  first  and  third  arguments  and 
non-decreasing  in  the  other  two.  Symbolically  J?(j,  T>  T)- 


We  now  give  some  examples  of  monotone  fluxes  H : 


(1)  Lax- Friedrichs  [2]: 


v+,v~)  =  H 


+  u 
2 


+  tr  \ 
2  ) 


~u~)-^ay(v+  -v  )  (2.3) 


where 


a1  =  max  \Hi(u,v)\, 

A<u<B 

C<v<D 


ay  =  max  \H2(u,v)\ 

A<u<B 

C<v<D 


(2.4) 


Here  Hi(u ,  u)  is  the  partial  derivative  of  H  with  respect  to  the  ith  argument.  The  flux 
HLF  is  monotone  for  A  <  u*  <  B,  C  <  v*  <  D; 


(2)  Godunov  type  [5]: 


Hg(u+,u~,v+,  v  )=  ext 

u€/(u~  ,u+) 


ext 

t»€/(v-  ,v+) 


HM 


(2.5) 


where 


I(a,b)  =  [min(a,6),  max(a,  6)] 


and  the  function  ext  is  defined  by 


ext  = 

«€/(a 


{min  if  a  < 

°<u<6 

max  if  a  > 


Ku<a 


(2.7) 


4 


As  pointed  out  in  [5],  since  in  general  minu  max„  H(u,  v )  ^  max„  minu  H(u,v),  we  will 
generally  obtain  different  versions  of  Godunov  type  fluxes  HG,  by  changing  the  order  of 
min  and  max.  Uniqueness  of  HG  happens  when,  e.g.,  H(u,v)  =  Hl(u)  +  H2(v ),  and  in 
many  other  cases.  Then,  by  [5],  — tHG(u,v )  is  the  exact  solution  to  the  Riemann  problem 
of  (2.1),  i.e.,  this  is  the  viscosity  solution  of  (2.1)  for 

(  u+, x  >  0  (  v* ,y>  0 

4>o(x,y)  =  xuo(x)  +  yuo(y),  uo(*)=<  _  U0(y)  =  <  _  (2.8) 

(  t£  ,x  <  0,  l  v  ,y  <  0 

evaluated  at  (i,  y)  =  (0,  0),  and  t  >  0. 

For  this  reason  all  monotone  fluxes  can  be  regarded  as  approximate  Riemann  solvers  in 
this  sense. 

(3)  Local  Lax-Friedrichs  (LLF)  [8]: 

ALLF(u+,u-,  v+,  V-)  =  £±^)-ia'(u+,  u-)(u+,  u~)~  Iq>+,  v~)(v+-,  r) 

(2.9) 

where 

a'*(tt+,u-)=  max  |/fi(u,u)|,  ay(v+,  u~)  =  max  \Hi{u,v)\  (2.10) 

«€/(«", «  +  )  v£l(v~  ,v+) 

C<v<D  A<n<B 

In  the  appendix  we  prove  that  HLLF  is  monotone  for  A  <  u±  <  B,  C  <  <  D;  Notice 

that  Hllf  has  smaller  dissipation  than  HLF . 

(4)  Roe  with  LLF  entopy  correction  [8]: 
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Hrf(u+,u  ,v+,v  )  =  ^ 


(  H(u*,v*)-,  if  Hi(u,  v)  and  H2{u,  u)  do  not  change  signs 
in  u  €  7(u~,u+),  v  €  I(v~,v+); 

H(u+$u  }v *)  —  |at(u+,u~)(u+  —  u~ );  otherwise  and  if  H2(u,  v) 
does  not  change  sign  in  A  <  u  <  B,  v  €  7(v~,  t>+) 
77(u*,  — )  —  |ay(u+,  u~)(i>+  —  u~);  otherwise  and  if  H\(u,  v) 
does  not  change  sign  in  u  €  7(u-,u+),C  <  v  <  7), 
l  HLLF(u+ ,u~  ,v+  ,v~);  otherwise 


(2-xl) 


where  ti*,  u*  are  defined  by  upwinding: 


u*  = 


u+  if  H\(u,  u)  <  0 
u~  if  H\(u,  v)  >  0 


v  = 


r>+  if  H2(il,  v)  <  0 
v~  if  H2(u,  v)  >  0 


(2.12) 


We  will  prove  in  the  appendix  that  HRF  is  monotone  for  A  <  u*  <  B,  C  <  <  D. 

Notice  that  ERF  is  simple  to  code,  purely  upwinding,  with  almost  as  small  dissipation 
as  HG .  Also  notice  that  HRF  is  not  continuous:  for  example  if  H(u)  —  (one  space 

dimension)  then  77FF(1,0)  =  0  but  HRF(  1,— e)  =  - j  •  1  •  (1  4-  e)  £ — —  |  ^  0. 

However,  this  type  of  discontinuity  does  not  hurt,  because  we  have 

|i7RF(tx+,  u~,  v+,  v~ )  —  Hg(u+,  u~,  <  M(\u+  —  u_|  +  |u+  —  r“|);  (2.13) 

hence  we  still  get  consistency  and  accuracy. 

REMARK  2.1.  A  flux  with  even  smaller  dissipation  than  HLLF  is 


H 


u 


,  V  ,  V 


-)  =  JJ(ld±«L  2l±2Ll)- lQV,y±)(«+-u"Hfts(ii±,r±)(r+-r-) 
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where 


az(u±,v±)=  max  jffi(u,v)j,  as'(u±,v±)=  max  j/f2(u,u)!.  (2.15) 

u£l(u~  ,u  +  )  u€/(u~,u  +  ) 

t>€l(v~  ,v  +  )  v€l(v~,v  +  ) 

Unfortunately  it  is  not  monotone:  for  example,  if  H(u,  v)  —  eu+t;,  then  HLLLF( 2, 0, 0, 20)  > 

ffLLLF(0,0,0,20).  For  separable  Hamiltonians  with  H(u,v)  —  H1( u)  +  H2(v)  we  have 
ftLLLF  _  jftLLF 

We  now  begin  the  description  of  our  high  order  ENO  schemes.  Monotone  fluxes  play 
the  role  of  “building  blocks”.  The  ENO  interpolation  idea  in  (1.3)  is  used  to  obtain  high 
order  non-oscillatory  approximations  to  u±  =  <f>f  and  =  <j>± .  These  values  are  then 
put  into  a  monotone  flux  H(u+ ,  u~,  t+,  v~).  Time  accuracy  is  obtained  by  a  class  of  TYD 
Runge-Kutta  type  time  discretizations  [7].  We  summarize  the  algorithm  as  follows: 

Algorithm  2.1: 

(1)  At  any  node  (i,  j),  fix  j  to  compute  along  the  i-direction,  by  using  (1.3),  obtaining 


d_  p<t>.r 
dz  rx±\,j 


(a-.)- 


(2.16) 


Similarly  for  vf~.  Then  let 


Lij  =  -A 


(2.17) 


(2)  obtain  from  4>n  by  the  following  Runge-Kutta  type  procedure: 

k- 1 

k  =  l,...,r  (2.18a) 


(2.18b) 


We  can  take  r  =  r  and  positive  ate  and  for  up  to  third  order  r  <  3.  The  method 


(2.18)  can  be  proven  TVD  under  the  CFL  condition: 


A  =  ~  <  Cr  A0 

Ax  “ 


„  it. .i _ r _ _ i _ : _ /o  i  *t\  rr» _ j _ il_  _ i:*  •  _ 


(2.19) 


A!fe-rithm  2.1  is  formally  uniformly  r-th  order  in  space  and  time  in  smooth  regions 
(measured  bv  local  truncation  errors). 

Notice  that  in  the  algorithm  above,  we  need  to  evaluate  two  polynomials  P±^}  to  get 
u i.  If  the  monotone  flux  is  purely  upwind  and  there  is  no  “sonic  point”  (points  at  which 
Hi  or  H 2  changes  sign),  one  of  u+  and  u~  is  never  used.  We  thus  recommend  the  following: 


Algorithm  2.2: 

(1)  Compute  u*  =  and  t>J-  =  ■  If  Hi(u,  t>)  and  H2{u ,  v)  do  not  change  signs 

in  u  e  I{u~,  v+j),  v  €  /(t>~,  t-^),  then  compute  only  u’j  and  v*}  by  (2.16)  where  u*,  t>*  are 
defined  by  (2.12);  and  take  Lt]  =  —AtH(u*,v*);  otherwise  take  (2.17). 

(2)  Same  as  step  (2)  in  Algorithm  2.1.  D 

Notice  that  Algorithm  2.2  is  NOT  equivalent  to  Algorithm  2.1  with  H  —  HRF .  Since 
we  expect  sonic  points  to  be  isolated,  Algorithm  2.2  is  usually  almost  twice  as  fast  as 
Algorithm  2.1. 


Rf.MaRK  2.2.  Notice  that,  in  smooth  regions,  by  Taylor  expansion, 


d_ 


A ±<i>,  _  Ax 
Ax  "  T 


If  we  choose,  instead  of  (2.16), 


where  the  projection  P  is  defined  by 


P[a,6](y)  =  { 


y,  if  a  —  b<y<a  +  b 
a  —  b,  if  y  <  a  —  b 

a  -f  6,  if  y  >  a  +  6 


(2.21) 


(2.22) 


(2.23) 
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we  will  still  have  uniform  high  order  accuracy  u*  =  {<pz  )XJ  -f  0(  Axr )  in  any  region  where 
l^nl  5:  2M.  Algorithm  2.1  will  then  give  a  scheme  which  deviates  from  a  monotone 
scheme  by  M  At  Ax,  hence  we  trivially  obtain  convergence  to  the  viscosity  solution  through 
the  theory  for  monotone  schemes.  In  practice  we  do  not  recommend  (2.22),  because  the 
parameter  M  is  not  intrinsic  -  it  has  to  be  adjusted  for  each  individual  problem.  See  [7, 
p.  452]  for  a  discussion  of  a  similar  situation  for  conservation  laws. 

REMARK  2.3.  When  implementing  (1.3)  we  use  undivided  differences: 

V*0',0)  =  ^  (2.24a) 

vU,k)  =  vU  +  l.k-l)-v(j,k-l),  k  =  l,...,r+l  (2.24b) 

The  computation  of  (2.24)  can  be  easily  vectorized.  The  ENO  stencil-choosing  process  is. 
for  computing  u+  =  (^x)+,  starting  with  i(j)  =  j  and  performing 


if  (abs [<p{i(j),k)).  gt-  abs(^(i(j)-l,*)))  i(j)  =  i{j)  -  1  (2.25) 

for  k  =  2, . . .  ,r,  where  i(j)  is  the  left-most  point  in  the  stencil  for  P’r,r,  (r).  This  can  also 

Jt  y 

be  vectorized.  Finally 


5^c(*0‘)  -j\k)?(i(j),k) 


.26) 


k=  1 


where 


.  m  +  fc  —  1  m-t-fc  — 1 

C(m,*)=-7  yi  n  (“*)• 


(2.27) 


*=m  (=m 


Notice  that  the  small  matrix  c  is  independent  of  y>,  is  only  computed  once  and  then  stored. 
(2.26)  can  be  vectorized  easily  as  well. 


3.  NUMERICAL  RESULTS 
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Example  1.  One  dimension.  We  solve 


i 

) 

l 

! 


i 


<f>t  4-  H(4>z)  —  0 
^(x,0)  =  —  cos  7rx 


-  1  <  X  <  1 


with  a  convex  H  (Burgers’  equation): 


H(u) 


(u  +  a)2 

2 


(3.1) 


(3.2) 


and  a  non-convex  H : 

H(u)  =  —  cos(u  +  a).  (3.3) 

Notice  that  if  we  let  v  =  <f>x  +  a,  f(v)  =  H(v  —  a),  then  (3.1)  becomes  a  conservation  law 

f  Vt  +  f{v)z  =  0 

<  -  1  <  x  <  1  (3.4) 

(  v(x,  0)  =  a  +  tt  sin  7rx 

which  is  a  standard  test  problem  for  conservation  laws  (e.g.  [7]).  We  can  easily  use  the 
method  of  characteristics  to  obtain  the  exact  solution  of  (3.1)  through  that  of  (3.4). 


We  take  a  —  1  and  compute  the  result  to  t  =  f  i  =  (when  the  solution  is  still  smooth) 
and  to  t  =  =  Ipr  (when  the  solution  has  a  discontinuous  derivative).  We  print  out  the 

L\  and  L0 c  errors,  in  Table  3.1,  for  selected  first  order  monotone  schemes  and  third  order 
ENO  schemes  in  smooth  regions,  i.e.  the  whole  region  [-1,1]  for  t  =  t\  and  the  region 
\x  —  x,|  >  0.1  for  t  =  f 2  where  x,  is  the  location  of  any  discontinuity  of  the  derivative.  We 
also  present  the  graphs  of  the  numerical  solutions  (in  diamonds)  versus  the  exact  solutions 
(in  solid  lines)  in  Figures  1  and  2. 
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Table  3.1  and  L»  Errors  in  Smooth  Regions  for  (3.1) 


-  0.5/tt2 


-  1.5/tt2 


Godunov 


.ETO-3-LF 


EH) -3 -God 


1 

10 

20 

1»0 

80 

10 

20 

4o 

i.oi(-i) 

5.13-2) 

2  59(-2) 

1.31(-2) 

2.01( -1) 

1.13(-1) 

6.13-2) 

3 

1.83/  -1) 

1.03(-1) 

6.03(-2) 

3.23(-2) 

3-17(-l) 

1.72(-1) 

8.43-2) 

3-33-2) 

1.50(-2) 

7.63-3) 

3-87(-3) 

5.93-2) 

3-13-2) 

l.?o(-2) 

7.43-2) 

3.6l{  -2) 

3&ES 

1.01(  -2) 

1.13-1) 

5  *93(  -2) 

3.00(-2) 

1.13(-2) 

1.83(-3) 

2.59(A) 

U.27( -5) 

1.65(-2) 

2.33  -3) 

2.63  A) 

3.13-2) 

5.08(-3) 

8.75(A) 

2.07(A) 

3-40(-2) 

7.8l(-3) 

1.33(-3) 

5.7H-3) 

9.55(A) 

1.41(A) 

2.53-5) 

7.37(-3) 

1.13-3) 

1.43  A) 

1.52(-2) 

4.43-3) 

7.55(A) 

1.53  A) 

2.29( -2) 

6.73-3) 

9.79(A) 

*1 

5.79(-3) 

9-53  A) 

1.43 A) 

2.53-5) 

8.19(-3) 

1.15(  -3) 

1.43  A) 

Q 

1-53-2) 

4.43-3) 

7.55(A) 

1.52(A) 

2.29( -2) 

5.93-3) 

9.97(A) 

H 

2.70(  -2) 

1.19( -2) 

5.53-3) 

2.53-3) 

5.5K-2) 

2.83(-2) 

1.33(-2) 

Is. 

4.i5(-a) 

2.33-2) 

1.33-2) 

6.93-3) 

1.03C-1) 

4.73-2) 

2.29(-2) 

B 

1.52(-2) 

5.89(-3) 

2.43-3) 

1.13 -3) 

2.4l( -2) 

1.13(-2) 

5.45(-3) 

a 

3-53-2) 

1.22(  -2) 

5.88(-3) 

2.57(-3) 

4.69(-2) 

2.53  -2) 

1.23-2) 

Ui 

1.63-3) 

5.63  a) 

8.93-5) 

9A7(-*) 

8.8l(-3) 

1.35(-3) 

1.84(A) 

S3 

6.23(-3) 

2.33( -3) 

4.81(A) 

1.03  A) 

1.39(-2) 

9-07(-3) 

1.53(  -3) 

B 

1.05(-3) 

3.17(A) 

5.93-5) 

7.57(-6) 

7.05(-3) 

9.10(A) 

1.51(A) 

m 

3.13(-3) 

1.17(-3) 

2.63  -10 

4.39(-5) 

1.3U-2) 

5*93(  -3) 

1.39(-3) 

B 

1.1K-3) 

3.21(A) 

5.93-5) 

7.57(-6) 

6.3K-3) 

9.10(A) 

1.53  A) 

S3 

3.WK-3) 

1.17(-3) 

2.63-4) 

*.39(-  5) 

1.27(-2) 

5.93-3) 

1.39(-3) 

80 


3 . 16(  -2  ) 


4  53-2) 


8.50(-3) 


l*5l(  ”2) 


4.l7(-5) 


1.07(A) 


2.05(-5) 


1.04(-4) 


2.85(-5) 


1.0U(-4) 


6.U3(-3) 


1.12(-2) 


2.60(-3) 


5.76  -3) 


2.U3(-M 


2.14(-L> 


2.53-5) 


2.0l(A) 


2.6o(  -5) 


2.01(A) 


Remark  3.1.  From  Table  3.1  and  Figures  1  and  2  we  can  observe  that: 


(i)  the  resolution  of  third  order  ENO  schemes  with  10  points  is  roughly  the  same  as  that 
of  the  corresponding  monotone  schemes  with  80  points; 

(ii)  ENO-3-Godunov  and  ENO-3-RF  have  roughly  the  same  resolution,  even  if  the  latter 

is  much  simpler  than  the  former  and  only  takes  about  half  time.  □ 


EXAMPLE  2.  Two  dimensions.  We  solve 

j  4>t  +  h(<j>z,<f>y)  =  0 

1  4>(x,y,  0)  =  -cosir^j^j 
with  a  convex  H  (Burgers’  equation): 


-2  <  x,y  <  2 


H(u,v ) 


(u  +  v  4-  a)2 

2 


(3.5) 


(3-6) 


and  a  non-con  vex  H : 

H(u,  v)  =  —  cos(u  +  v  4-  a).  (3-7) 

Notice  that,  under  the  transformation  f  q  =  (3.5)  -  (3.6)  -  (3.7)  become  (3.1) 

-  (3.2)  -  (3.3)  in  the  £  direction.  We  can  thus  use  the  one  dimensional  exact  solution  to 
analyze  our  numerical  results.  Since  we  use  (x,  y)  coordinates,  this  is  a  true  two  dimensional 
test  problem. 

We  again  take  a  =  1  and  compute  to  t  =  ti  =  ^  and  t  =  t2  —  ^r-  Some  results  are 
presented  in  Table  3.2  and  Figure  3. 
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■-eoa(  u+r+lf^ 


Table  3.2  L\  and  Loo  Errors  in  Smooth  Regions  for  (3.5) 


t  »  0.5/n 


t  «  1-5/tt 


ETO-3-Ur 


1 

@ 

l» 

8o8 

102 

mm* 

■■ 

1.01( -1) 

5.10<-2) 

2.59C-2) 

1.3K-2) 

2.24(-l) 

1.25(-1) 

6.79(  -2) 

1.8K-1) 

1.08(-1) 

6.03(-2) 

3.23(-2) 

3 .17(  -1) 

1.72(-1) 

8.98(-2) 

1.15(*«) 

1.83(-3) 

^.59(-4) 

4.25( -5) 

1.83(-2) 

2.58(-3) 

2.90(-4) 

3.37(-S) 

5.09(-3) 

0-75( -4) 

2.1K-4) 

3.4l(-2) 

7.82(-3) 

1.34(-3) 

2.70(-S) 

1.19(-2) 

5.56< -3) 

2.58(-3) 

8.01(  -2) 

3-56(-2) 

1.67( -2) 

4.1 5(-«) 

2. 33-2) 

1.30(-2) 

6.93  -3) 

1.03( -1) 

4.72( -2) 

2.29(-2) 

1.65(-3) 

5.68(-4) 

8.93(-5) 

9-23( -6) 

1.03( -2) 

1.69(-3) 

2.20)-4) 

6.23(  -3) 

2.33(-3) 

4.82(-4) 

1.07( -4) 

1.42( -2) 

9-10(  -3) 

1.56(-3) 

8o2 


3-51(-2) 


4.59(>2) 


4.6U(-5) 


1.15(  -4) 


8.03(-3) 


1 . 12(  -2) 


4.93(  -5) 


6.13(-4) 


Remark  3.2 

(i)  By  comparing  Table  3.2  with  Table  3.1  we  can  see  that  ENO  schemes  perform 
equally  well  in  two  dimensions; 

(ii)  Notice  that,  except  for  a  sharper  discontinuity-in-derivative  resolution,  we  cannot 
see  much  difference  between  Figures  3(a),  3(c)  (first  order  monotone  schemes)  and 
3(b),  3(d)  (third  order  ENO  schemes).  However,  from  Table  3.2  we  can  clearly  see 
a  large  difference  in  the  resolution  of  the  solution  in  smooth  regions.  This  indicates 
the  limitations  of  graphical  presentations; 
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(iii)  In  this  two  dimensional  case,  the  Godunov  flux  is  considerably  more  complicated 
to  program  than  LF  or  RF ,  with  a  not-so-significant  improvement  in  resolution  for 


ENO-3. 

EXAMPLE  3.  We  solve  a  two  dimensional  non-convex  Riemann  problem 

{4>t  +  sin(<j>x  +  4>y)  =  0 
<f>(x,y,0)  =  7f(jy|  -  |x|) 


□ 


(3.8) 


to  investigate  the  resolutions  of  different  building  blocks,  the  behavior  of  different  versions 
of  Godunov  flux  (2.5);  and  convergence  to  viscosity  solutions.  The  results  are  in  Figure  4 
and  5:  From  the  graphs  and  computer  outputs  we  can  observe: 

(i)  ENO-3  with  Gl,  G2  (two  versions  of  Godunov  fluxes),  LF  and  RF  as  building 
blocks  are  all  convergent  to  the  viscosity  solution,  with  a  much  sharper  resolution 
for  the  discontinuities-in-derivative  than  the  corresponding  first  order  monotone 
schemes; 

(ii)  ENO-3-RF  has  roughly  the  same  resolution  as  EN0-3-Godunov,  with  a  much  sim¬ 
pler  program  and  a  reduced  computer  cost; 

(iii)  The  difference  between  two  versions  of  Godunov  fluxes  is  very  small:  the  average 
difference  at  t  =  1  is  around  1000  times  smaller  than  the  L\  errors. 

(iv)  ENO-3,  using  the  Roe  flux  as  a  building  block  without  entropy  corrections,  i.e., 
Algorithm  2.2  without  using  (2.17)  for  entropy  corrections  in  “sonic  cells”,  converges 
to  an  incorrect  solution  just  as  the  first  order  Roe  scheme  (Figure  5).  This  indicates 
the  importance  of  entropy  corrections  in  “sonic  cells”. 

EXAMPLE  4.  We  solve  the  following  problem  related  to  control-optimal  cost  determi¬ 
nation: 


{4>t  -  (sin  y)(f>z  4-  (sin  x  -I-  sign(  4>y  ))4>v  -  \  sin2  y  -  ( 1  -  cos  x)  =  0 

<t>(x,y,  0)  =  0 


(3.9) 
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assuming  periodicity.  The  results  at  t  =  1  are  presented  in  Figure  6.  Notice  that  third 
order  ENO  schemes  have  sharper  discontinuities-in-derivative  resolution  than  first  order 
monotone  schemes.  For  this  problem  the  interesting  quantity  is  the  optimal  solution  w  = 
sign(<£v),  Figure  5d.  A  sharper  discontinuities-in-derivative  resolution  means  smaller  error 
for  w  in  the  neighborhood  of  each  such  point. 


4.  CONCLUDING  REMARKS 

We  have  generalized  ENO  schemes  for  conservation  laws  to  Hamilton- Jacobi  type  equa¬ 
tions.  Computational  results  indicate  good  accuracy  in  regions  of  smoothness,  sharp  dis¬ 
continuities  in  derivatives,  and  convergence  to  the  correct  viscosity  solutions.  Algorithm 
2.2  (ENO-RF)  is  usually  preferred,  due  to  its  simplicity  to  program,  reduced  computa¬ 
tional  cost,  and  its  excellent  resolution,  which  is  comparable  to  the  results  using  the  much 
more  complicated  Godunov  type  building  blocks. 
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APPENDIX 


i 

We  prove  that  HLLF  and  HRF  are  both  monotone.  To  simplify  the  exposition  we  only 
consider  the  one  dimensional  case.  The  proof  for  the  multi  dimensional  case  is  similar.  i 

LEMMA  A.l  HLLF  is  monotone. 

Proof.  Hllf  in  one  dimension  is  defined  by 


HLLF(u+,u-)  =  H 


(^) 


(A.l) 


We  assume  it*  >  and  want  to  prove  HLLF (uf  ,u  )  <  HLLF (u* ,u  ).  Let  D  = 
Hllf(u* ,  u~)  —  HLLF(u2  ,  u~).  This  equals: 


H 


-  max  \H'{u)\(u+-u  )+-  max 

^  u€/(u~  ,uf)  v€I(v~  ,uf) 


CASE  (i)  >  u  .  We  have  for  u  <  ■  <  uf ; 


-Uj)-  max  |jy'(tt)|(«+ -  tt  )+  max  \H'(u)\(u2 

u~ <u<u+  u_  + 

<  r|#'(0(u*  ~  u2  )  -  max  !H'(u)|(u^  -  u“)  +  max  )if'(u)|(u^ 

^  L  u~  <u<ujt  u-  <u<u  + 


=  ”  u2  ) 


H'(0~  max  |ff'(u)| 


u ~  £u<uj*’ 


<0; 


Case  (ii)  u  >  u+  >  ,  similar  to  case  (i); 


-u  ) 
-  u~) 
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CASE  (iii)  uF  >  u  >  u^.  We  have,  for  it  <  £  <  ■>~^u  -  <  it4  and  it  4  <  <  t]  <  u  , 


D  =  \  H'(0(ut  -  u  )  +  H'(rl )(u  -ut)~  max  I H'(u)\(uF-u  ) 

2  U-<U<li  + 


+  max  |iif '(tx)|(tx4  —  u  ) 


U2  —U—U~ 


=  l;(ut  ~  u  )  H'(0~  max  )H'(ii)j  +  ^-(u  -u^)H\r])-  max  |i/'(u)|  <0 

2  U_<U<U*  2  tij<U<U~ 


Hence  we  proved  H(|,  •).  We  can  similarly  prove 


Lemma  A. 2  HRF  is  monotone. 


Proof.  HRF  in  one  dimension  is  defined  by 


H(tx4)  if  H'{u)<Q  in  u  £  I{  u  ,tx4) 


HRl  (ix4,u~)  =  l  H(u~)  if  H'(u)  >  0  in  u€J(u-,u+)  (A.2) 

l  HLLF(u+  ,u~)  otherwise. 

We  assume  uf  >  uf  and  want  to  prove  HRF(uF ,  u^)  >  HRF( tx4,  tx^"). 


Let  D  =  HRF(u+,uf)  -  Hrf(u+ ,  tx 2“). 

CASE  (i)  HRF(u+ ,  uf )  =  tfG(tx4,tx-)  and  RF(tx4,  tx")  =  JJG(tx4,  u~ )  or  HRF{u+,  u~ ) 
Hllf(uf  ,  txj~)  and  HRF(u+,  uf)  =  HLLF(u+ ,  tx^),  then  automatically  D  >  0. 


Case  (ii)  HRF(tx+,uD  =  tf(u-),  JfRF(tx4,uJ)  =  tfiLF(tx4,  txj  ).  Then  tf'(tx)  >  0  in 
J(uf,u+)  but  H'(u )  changes  sign  in  J(tx^,tx4),  hence  J(tx^,tx4)  £  J(tx|~,  tx4).  Therefore 
tx4-  >  txr .  We  then  have  either 
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(a)  u+  >  uf  >  uj,  and,  for  u2  <  (  <  u+  we  have 
D=\  H'(0( 2u~-u2-u+)+  max  \H'(u)\(u+  -  u2) 

2  _  U~ <tl<U  + 

=  i(u+-u~)  max  \H'(u)\  -  H'(£)  max  |ff'(u)|  +  H'(Q  >0 

2  L«J<«<“+  J  2  Lu~<u<u+  J 

or  (b)  uf  >  u+  >  uj,  and,  since  H(u^)  >  H(u+ )  due  to  the  fact  that  H'(u )  >  0  in 
[u+,  uj“],  we  have,  for  u+  >  £  >  -  ^Uz-  >  uj" , 

D  =  ff(uf)  -  J?iLF(u+,t^)  >  H(u+)  -  ^iLF(u+,u-) 

=  ^  )  +  max .  i-ff,(u)i(u+  -  u2 ) 

2  <ti<Uj 

=  I(u+_u-)f  max  |ff'(u)|  +  H'(0)  >  0; 

2  \  Uj  <u<u+  / 

Case  (hi)  ifRF(u+,u;-)  =  H(u+),  HRF(u+,u2)  =  HLLF(u+,u~) 

As  in  case  (ii)  we  can  again  deduce  u+  >  u2',  hence,  for  u2  <  ^  a  <  £  <  u+,  we  have 

D  =  ff(u+)-^LiF(u+,«“) 

=  i  H'(£)(u+  -  u2)  -  max  |jy'(u)|(u+  -  u2)  >0; 

2  _  Uj<U<U+ 

Case  (iv)  HRF(u+,  uf )  =  tfLLF(u+,  uf ),  HRF(u+,u2  )  =  (u^ )  similar  to  case  (ii); 

Case  (v)  Hrf(u+,u~ )  =  tfLiF(it+,uf ),  tfRF(u+,uJ)  =  tf(u+)  similar  to  case  (iii). 

We  have  proved  #(•,  T)-  Similarly  for  H(l,  •)  O 
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lc  :  First  order  Lax-Fnedricbs,  80  points 


-0.5  0.0  0.5 


le  :  Third  order  ENO-Godunov,  10  points 


Id  :  Third  order  ENO-LF,  10  points 


If  :  Third  order  ENO-RF,  10  points 


Figure  1  :  One  dimensional  Burgers’  equation 
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3c  :  First  order  Lax-Friedrichs,  non-convex  H 


3d  :  ENO-3-LF,  non-convex  H 


Figure  3:  Two  dimensions,  40  x  40  points 


4a  :  First  order  Lax -Friedrichs 


4b  :  ENO-3-LF 


4c  :  First  order,  Godunov,  version  I 


4d  :  ENO-3-G1 


Figure  4:  Riemann  problem  (3.8),  40  x  40  points,  t  =  1. 
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